Vortex formation in a rotating two-component Fermi gas 
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A two-component Fermi gas with attractive s-wave interactions forms a superfluid at low tem- 
peratures. When this gas is confined in a rotating trap, fermions can unpair at the edges of the 
gas and vortices can arise beyond certain critical rotation frequencies. We compute these critical 
rotation frequencies and construct the phase diagram in the plane of scattering length and rota- 
tion frequency for different total number of particles. We work at zero temperature and consider a 
cylindrically symmetric harmonic trapping potential. The calculations are performed in the Hartree- 
Fock-Bogoliubov approximation which implies that our results are quantitatively reliable for weak 
interactions. 



I. INTRODUCTION 

A characteristic feature of superfluids is the appear- 
ance of vortices when they are rotated. This fact has 
been used to demonstrate that a two-component atomic 
Fermi gas becomes a superfluid at sufficiently low tem- 
peratures . A superfluid state can be created in such 
a gas by trapping fermionic atoms in two distinct hyper- 
fine states. The interaction strength between the two 
components can be controlled by an external magnetic 
field. When the interactions are tuned to be attractive, 
and the atoms are cooled to sufficiently low temperatures, 
the components will form pairs via the Cooper instabil- 
ity. Due to this pair formation the Fermi gas becomes 
a Bardeen-Cooper-Schrieffer (BCS) superfluid as was en- 
visaged in Refs. H, Q. 

The response of the superfluid to rotation can be inves- 
tigated by rotating the trapping potential with a certain 
frequency fl. For a non-rotating trap the entire gas will 
form a superfluid without vortices. Let us now imagine 
increasing the rotation frequency at zero temperature. 
Up to a certain critical rotation frequency, the superfluid 
will stay in the vortex-free state carrying zero angular 
momentum. For low temperatures the angular momen- 
tum will be quenched, as has also been observed exper- 
imentally Above the critical frequency, angular mo- 
mentum will be inserted in the gas by either unpairing 
the fermions near the edges of the gas, by formation of 
vortices, or by the combination of both effects. The goal 
of this paper is to compute the rotation frequencies at 
which these transitions take place. 

Besides the experimental investigations [H-ll, 0| : vari- 
ous theoretical studies of rotating two-component Fermi 
gases have been performed (for a review see e.g. Ref. [1]). 
The profile of a single vortex was analyzed in this con- 
text for the first time in Ref. [Qi] using a Ginzburg-Landau 
approach and in Ref. by solving the Bogoliubov-de 
Gennes equation. In Ref. it was concluded that a 
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single vortex can induce a sizable density depletion at its 
core. This was also found to be the case for a vortex lat- 
tice [l2j . Density depletion is important, since it allows 
the vortices to be detected experimentally [l| . The vortex 
profile was investigated in a population imbalanced gas 
in Ref. [l3j and in situation in which the two components 
have unequal mass in Ref. [13] • Real-time dynamics of 
vortices has been studied in Ref. [Tsj . 

Vortex lattices in two-component Fermi gases have 
been examined in several ways in Refs. [l^, [3 [13 ■ At 
high rotation frequencies a completely unpaired phase is 
preferred over a vortex lattice. The critical rotation fre- 
quency corresponding to this transition was computed in 
Refs. [H, [l^. A vortex lattice can also be destroyed by 
heating the gas. The corresponding critical temperature 
was computed in Refs. [H, [13] ■ When the number of 
trapped components is unequal, a vortex lattice might 
be formed within the Fulde-Ferrell-Larkin-Ovchinnikov 
(FFLO) phase. Its melting temperature was investigated 
in Ref. [2l| and its upper critical rotation frequency in 
Ref. [12. 

The first analysis of the critical frequency fic for vortex 
formation in a two-component Fermi gas was performed 
in Ref. |23|]. To obtain ilc the Helmholtz free energy dif- 
ference AF between a vortex with unit angular momen- 
tum located at the center of the trap and the vortex-free 
superfiuid was estimated at 17 = 0. For a cylindrically 
symmetric infinite well as the trapping potential, AF was 
obtained by solving the Bogolibuov-de Gennes equation 
in Refs. [10, 24]. The free energy difference arises from 
the loss of condensation energy at the vortex core, the 
kinetic energy of fermions circulating around the vortex 
core, and the energy needed to expand the cloud to ac- 
commodate the excess particles removed from the vortex 
core (the latter effect was not considered in Ref. [1^ , but 
was taken into account in Refs. [i3;[13)- When rotating 
the trap, the free energy decreases with flLz, where Lz 
is the angular momentum contained in the gas. For the 
vortex- free superfluid Lz ~ 0, while for the vortex with 
unit angular momentum Lz — Nh/2, with N the num- 
ber of trapped particles. These estimates are correct if 
rotating the trap does not cause unpairing near the edges 
of the gas. The critical rotation frequency in this case is 
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flc = 2AF/{NH). This frequency is similar to the lower 
critical magnetic field in type-II superconductors p5| . 

Fermions confined in a rotating trap can unpair near 
the edges |2 6|-[2 8l|. This effect was not considered in 
Refs. [lOl I23I, [2JI . In this paper we will take into account 
this possibility in order to obtain a more reliable value of 
ilc- We will make detailed study of how ilc varies with 
the interaction strength and the number of trapped par- 
ticles. Furthermore we will compute the critical rotation 
frequency for unpairing. We consider a cylindrically sym- 
metric harmonic trapping potential. Our calculations are 
carried out in the Hartree-Fock-Bogoliubov approxima- 
tion. Therefore, we expect that our results are reliable 
only in the weak coupling limit. 

Let us finally remark that vortices have also been ob- 
served in rotating Bose-Einstein condensates (EEC) of 
bosonic atoms [29| and of bosonic dimers composed of 
fermionic atoms (for a review see e.g. Ref. j3Q| ) . The 
behavior of vortices through the BEC-BCS crossover was 
investigated experimentally in Ref. [H . In this article we 
only discuss the BCS regime. A computation of the crit- 
ical rotation frequency for vortex formation in a BEG is 
discussed in Ref. [3l|. Theoretical studies of behavior 
of vortices through the BEC-BCS crossover have been 
performed in Refs. [32| - [33 |. 

This article is organized as follows. In Sec. |ll]we will 
introduce the action from which one can derive the prop- 
erties of the Fermi gas. To achieve this in practice we will 
adopt the two-particle irreducible (2PI) effective action, 
which we explain in Sec. IIIIl From the 2PI effective ac- 
tion one obtains the Dyson- Schwinger equation which is 
the main equation we have to solve. This can be achieved 
by finding the solution of the Bogoliubov-de Gennes 
equation, which is explained in Sec. |IVl The numeri- 
cal methods by which we have solved the Bogoliubov-de- 
Gcnnes and the Dyson-Schwinger equation are discussed 
in Sees. Ivl and I VII respective Iv. The reader who is not in- 
terested in the details of the calculation can immediately 
go to Sec. I VIII where we present the results. The phase di- 
agrams presented in Figs. I14land ll5l are our main results. 
We draw our conclusions in Sec lVIIll Several details are 
relegated to the appendices. In Appendix we review 
the 2PI effective action. A derivation of the Bogoliubov- 
de Gennes equation is presented in Appendix [BJ To solve 
the Bogoliubov-de Gennes equation numerically, we will 
use a basis based on Maxwell polynomials. We discuss 
the computation of the quadrature weights and nodes of 
these polynomials in Appendix[Cj In Appendix [P] we de- 
rive the representation of the single particle Hamiltonian 
in the basis we will employ. 



II. SETUP 



idity is driven by attractive s-wave interactions. The in- 
teractions among the same components can be neglected 
since the Pauli principle admits s-wave interactions only 
between the different species. We will denote the s-wave 
interaction potential as V{x) and specify it below. Under 
these assumptions the interacting Fermi gas is described 
by the following action (see e.g. Ref. (s^), S = Skin + S'int 
where 
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a=t,i (2) 

X V{x - y)d{T, - Ty)^MY)MX). 

Here ipaiX) is the (path-integral) quantum field corre- 
sponding to the a component, and X = {x,Tx). We 
write integration over spatial coordinates and imaginary 
time T as J dX = J^^ dr^ J d^x with the inverse temper- 
ature /3 = l/{kBT). Here and in the rest of the article 
X± = {x,Tx ±7]) and 77 is an infinitesimal small posi- 
tive number. We have inserted X+ and 1+ in Eq. ^ 
in order to maintain the correct ordering of the fields 
in the path-integral. We will achieve this in a different 
way for the kinetic term and explain this at the end of 
Appendix |B] Furthermore — f is equivalent to J,, /i^ 
denotes the chemical potential, and H{il) is the single- 
particle Hamiltonian. We will assume that the particles 
are trapped in a potential that is rotating in the x-y 
plane with angular frequency fl. It is then convenient 
to perform the calculation in the rotating frame. The 
single-particle Hamiltonian H(Cl) in the rotating frame 
reads (see e.g. Ref. jssj ) 
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+ U{x) - flL^ 



(3) 



where M is the fcrmion mass and = xpy — ypx is the 
z-component of the angular momentum. The trapping 
potential U (x) realized in experiments is typically har- 
monic. In this paper we will study a cylindrically shaped 
trap given by the potential 



U{x) = ^Muj^{x^ + y^), 



(4) 



which implies confinement in the x-y plane, and infinite 
extension in the z-direction. In an experiment this regime 
can be reached by choosing the trapping frequency in 
the z-direction much smaller than in the x-y-direction. 
In cylindrical coordinates, x = (p cos 0, p sine/), z), the 
single-particle Hamiltonian reads 



Let us consider a two-component Fermi gas in which s- 
wave interactions are dominant and label its components 
by a =t, i- Typically in experimental realizations the 
higher partial waves can be neglected and the superflu- 
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where = -~ihd/d4>. The normahzed eigenfunctions 
ipnmp^ (x) of this Hamiltonian can be written as a product 
of three functions, 



{P) fruity' 



(6) 



with the radial quantum number n — 0, 1, . . ., the an- 
gular momentum quantum number m G Z, and the mo- 
mentum in the z-direction pz = 2Trhnz/L with Uz G Z. 
The constant L denotes the length of the system in the 
z-direction. In this article we will consider the limit 
L cx). Then x En, = Jdpz/{2Trh). The function 
fm{<l>) is a normalized eigenfunction of Lz and is given 
explicitly by 



/n^(0) 



1 



2tt 



(7) 



The radial eigenfunctions are given by 

where L'^{x) denotes the generalized Laguerre polyno- 
mial which has degree n. Furthermore p = p/X with the 
harmonic oscillator length A = (fi/Mw)^/^. Normaliza- 
tion gives c^^^j — 2n!/[A^(n-|-|rn|)!]. The energy spectrum 
of H(n) is given by 



hoj{l + 2n+\m\) ~ hnm+ j^. (9) 



At low enough temperatures and densities, the typi- 
cal wavelength of the particles will be much longer than 
the range of the interaction potential. In that case the 
detailed structure of the potential is unimportant and 
the only relevant interaction parameter is the scattering 
length a. To perform a calculation in this situation, one 
can just choose the most convenient potential that has 
scattering length a. Following Ref. fSg] we will use the 
Huang- Yang potential [37 1 



V{r) 



g<5(r)^7-, 



(10) 



where the coupling constant g = ATrafi^ /M . Pairing be- 
tween fermions requires attractive interactions, that is 
a < and hence g < 0. Note that the Huang- Yang 
potential is not equivalent to an ordinary 5-function po- 
tential gS{r), because the derivative operator also acts 
on the fields in Eq. The advantage of the Huang- 
Yang potential is that all relevant physical quantities be- 
come automatically convergent [36| . Furthermore, the 
coupling constant does not have to be renormalized so 
that the foregoing relation between the scattering length 
and the coupling constant always holds [1^. In the case 
of an ordinary (5-function potential one will encounter di- 
vergences. This will require a regularization prescription 
and a renormalization of the coupling constant. 



In order to perform calculations it is convenient to 
rewrite the action in the Nambu-Gor'kov basis. For that 
purpose we introduce the Nambu-Gor'kov fields 

*w = (*:<5), ,11) 

and write the kinetic part of the action, Eq. ([T]) , as 
^kin ^-hJdX I dX'-^^X)Go\X,X')^{X'), (12) 

where the bare inverse Nambu-Gor'kov propagator reads 
Go\X,X') = 

_}_( h-l-^+H{n)~p^ 
h\ h^-H{nY + p^ 

x5{X-X'). (13) 

Here we used the fact that H{~Vl) — H{il)*. In the 
Nambu-Gor'kov basis the interaction part of the action 
becomes 

5int = - y" / dXdYdX'dY' ■9^X')aa'i'iX) 

t^J (14) 

X ^\Y')a^^-^{Y)V^{X,Y-X\Y'), 

where (7+ — diag(l,0) and a-_ = diag(0, 1). The poten- 
tial Vq is given by 

V^{X,Y-X\Y') = ^[il + a)S{T,, ^Ty)Vix' -y) 

+ {l^a)SiT,~Ty,)Vix - y')\s{X^~X')d{Y.^-Y'). 

(15) 

Because the Huang- Yang potential contains a derivative 
operator, we had to introduce several (5-functions in or- 
der to separate the potential operator from the quantum 
fields. 



III. 2PI EFFECTIVE ACTION 

To study the interacting Fermi gas we will compute the 
resummed propagator Gij{X,X') ^ -(^^(A)*]^^')) 
and the grand potential $g by using the two-particle 
irreducible (2PI) effective action 38|. Readers who are 
not interested in the details of this formalism can im- 
mediately continue with Sec. IIVI where we discuss the 
Bogoliubov-de Gennes equation which follows from the 
2PI effective action. 

The 2PI effective action is also known as the Cornwall- 
Jackiw-Tomboulis formalism and is equivalent to the 
Luttinger-Ward functional approach [3£|. It has been 
applied previously to inve stig ate pairing in atomic gases 
[40| and in quark matter |4l| . 

The main advantage of the 2PI effective action is that 
it is a functional method which generates the resummed 
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FIG. 1: Feynman rules for the propagator and the vertex. 
When connecting propagators to the vertex, the arrows of 
the propagators have to point in the same direction as the 
arrows in the vertex. For each closed fermion loop one should 
include a factor —1. 





FIG. 2: All 2PI diagrams contributing to r2[G] at order g 
and order g^. Note the direction of the arrows. 



Nambu-Gor'kov propagator and the corresponding grand 
potential. Here, we will truncate the 2PI effective ac- 
tion at order g, which is equivalent to the Hartree-Fock- 
Bogoliubov approximation. This leads to the well known 
Bogoliubov-de Gennes equation. Another advantage of 
the 2PI method is that any truncation can be systemat- 
ically improved straightforwardly by taking into account 
higher order diagrams. This is necessary when extending 
our results to a strongly coupled Fermi gas. 

The 2PI effective action reads (see Appendix |X] for 
details) ^ 

r[G] = -TrlogG-i-Tr(GoiG-l)+r2[G], (16) 

where r2[G] is the sum of all 2PI diagrams generated 
from ^int with propagators G. The interaction vertex can 
be directly read off from Eq. (ITU) . We have displayed the 
Feynman rules for computing the 2PI effective action in 
Fig. [H All diagrams contributing to r2[G] up to order 
g'^ are displayed in Fig. [2] 

By minimizing r[G] with respect to G one obtains the 
Dyson-Schwinger equation 



G — Gq Xj[G], 



(17) 



where the IPI self-energy is S[G] = ST2[G]/dG. The 
grand potential $g — log Z is equal to the minimum 
value of F [G] //3 and reads 

<i>G = -^TrlogG-i - ^Tr(E[G]G) + ^T^IG], (18) 

where G is now the solution of the Dyson-Schwinger 
equation, Eq. pT|) . 



In order to perform a practical calculation, one has to 
truncate r2[G] at some order. In this paper we will take 
into account the contributions of order g to r2 [G] , which 
are represented by the first two diagrams in Fig. [51 Such 
truncation leads to the Hartree-Fock-Bogoliubov approx- 
imation. This approximation is expected to give accurate 
results for small coupling g. In general the Hartree-Fock- 
Bogoliubov approximation can also be applied to strongly 
correlated systems if the interaction kernel is renormal- 
ized appropriately. Such renormalization entails a re- 
summation of ladder diagrams in the Lippman-Schwinger 
equation. Since we are interested in the weakly interact- 
ing BCS limit there is no need to do this. 

Applying the Feynman rules to the first two diagrams 
in Fig. [2] we find that to order g, r2[G] is given by 

^2[G] = -^ Yj [ dXdYdX'dY' Va,{X,Y;X' ,Y') 

' a = ± 

X tr [G{X, X')<7a] tr [G(r, Y')a^a] 
+ j^Yl [ dXdYdX'dY'Va{X,Y;X\Y') 



i^[G{X,Y'WG{Y,X')a-^] 



(19) 



The relative minus sign arises because the first diagram 
in Fig. [2] contains two closed loops whereas the second 
diagram contains only one closed loop when following 
the arrows. 

The self-energy S [G] can be obtained by differentiating 
Eq. (fT9|) with respect to G. This yields 



nG]{X,X') = ^ fTa / dydr'tr [G{YX)^-A 

a— ± 



xVc.{X',Y-X,Y') 

f dYdY' a^G{YXy-o. 

a=±-' 

X Va{X',Y;Y',X). 



(20) 



The first term is the Hartree self-energy, the second one 
is the pairing contribution. 

The last two equations can be simplified by comput- 
ing the traces and inserting the explicit expressions for 
V and V{r). One has then to act the Huang- Yang po- 
tential on the propagators. It can be shown [36| that the 
diagonal components of G(X, X') are finite in the limit 
r — where r = |a; — a;'|. For these components the 
Huang- Yang potential acts as an ordinary 5-function po- 
tential. We will use this fact to simplify the expressions 
involving the diagonal components of G. On the other 
hand the off-diagonal components of G will have in gen- 
eral a singularity proportional to 1 /r when r — >■ (36j . In 
such situations the full form of the Huang- Yang potential 
needs to be taken into account. 

For later convenience we will now define the pairing 
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field as 



A(a;) = j A\'V{x-x')Gn{x,T-x',T). (21) 

Following Ref. [sB] we will split the off-diagonal compo- 
nent of the Nambu-Gor'kov propagator into a singular 
and regular part 



C 



lim Gti(a;, r; a;', r) = - + G^f (X, X), 



(22) 



where C is some constant and the superscript reg indi- 
cates the regular part of G. By inserting the Huang- Yang 
potential in Eq. ([2T|) one can see that the pairing field 
does not contain any singularity 36] 



^{x)=gG\J{X,X). 



(23) 



The diagonal components of the propagator can be ex- 
pressed in terms of number densities. To find the relation 
between the propagator and the densities we use the def- 
inition of the propagator in terms of field operators 

G,j{X,X') = -{T^^,{X)i>]{X')) 

= e{T' -r){¥^{X')i>,{X)) (24) 

-e{T-T'){i>,{x)^{x')), 

here denotes time ordering in imaginary time. Hence 
the number densities of the two species are related to G 
as 

n^{x) = {^{x,T)^:^{x,T)) ^G^^{X,X+), (25) 
n^{x) = {^?^{x,t)^^{x,t))=-G^^(X,X^). (26) 

From Eq. ^ it follows that G'^^{X,X') = 
(^;(X')^t(^)) and G^^{X,X') = {^{X')^l{X)) . This 
can be used to show that G)^^;{X,X') = Gn{X',X)* 
which implies that 



A(a;)* = j A^x' V{x - x')Gi^{x,t-x\t) 



gG\'HX,xy 



(27) 



Here we used the fact that the singular part of 
G^l{X' ^X)* is a function of \x — x'\, so that x and x' 
could be interchanged without changing the result. 

We can now use the above definitions to simplify 
Eqs. (IT^ and (1^01) . After computing the traces we find 
that to order g, r2[G] is given by 



^T2[G]^g j A'xn^{x)n^{x) 



(28) 



+ / A^xGu{X,xy^{x). 



The IPI self-energy to order g is 
^[^K^,^') - \ ( '^^'^^ ) KX-X'). (29) 



It follows that Tr(E[G]G) = 2r2[G] -f 0{g^) so that the 
grand potential $g, Eq. (|18p . to order g becomes 



$G = --^TrlogG ^-g / d^xn^{x)n_i{x) 



J d''xGu{X,xrA{x) 



(30) 



Since the particle number in the trap is fixed we solve 
the equations for the chemical potentials to obtain the de- 
sired number of fermions in each hyperfine state. The ap- 
propriate thermodynamic potential is then the Helmholtz 
free energy, which reads 



F = <i>G + A^t^t + fJ-iXi, 



(31) 



where N^,i denote the total number of particles of a par- 
ticular species. Since we consider a cylindrically shaped 
trap, F and Nf,^ are proportional to the length of the 
trap in the z-direction L. Since L is taken to be infinite 
is convenient to consider instead the free energy and par- 
ticle number per unit of the harmonic oscillator length A 
in the 2;-direction. For this reason we define 



F 



t4 



L/X' 



(32) 



Furthermore we will define Af to be the total number 
of particles per unit of length in the z-direction, Af = 

Aft + Mi- 



IV. BOGOLIUBOV-DE GENNES EQUATION 

To proceed, we will insert the explicit expression for 
the IPI self-energy, Eq. into Eq. (H?]). This yields 

the Dyson-Schwinger equation for the Nambu-Gor'kov 
propagator, 

G-\X,X') = -^(^h^+'H}jS{X ^X'), (33) 

with 



n 



H{n)-fi-i+gni{x) 
A*{x) 



A{x) 

-H{n)*+^l-gnt{x) 



(34) 



To solve the Dyson-Schwinger equation, one first inverts 
both the left- and right-hand sides of Eq. (p3| . As ex- 
plained in detail in Appendix [B] this can be achieved by 
solving the Bogoliubov-de Gennes equation (i^ 



n 



Ui{x) 
Vi{x) 



E, 



Ui{x) 
Vi{x) 



(35) 



The functions Ui(x) and Vi{x) have to be normalized as 
J d^x [\ui{x)\'^ + \vi{x)\'^] = 1. Using the exphcit ex- 
pression of G, Eq. (IB4I) . and Eqs. (PS)) . (1^51) one can now 
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read off the densities 

n^ix) = Y,fiE^)\u^{x)\^ (36) 

i 

n^{x) = (37) 



where f{E) = [exp{/3E) + 1]"^ is the Fermi-Dirac distri- 
bution function. 

As follows from Eqs. and (|B4)) to obtain A(a;) we 
need to extract the regular part of the propagator 

Gnix,T;x',T)=Y,fiE^)u,{x)v*ix'), (38) 



in the limit x' x. To do so we will use the method 
proposed in Ref. [36| with the improvements suggested 
in Refs. gUil. 

The sum over all modes in Eq. (|55| is logarithmically 
divergent for x = x' . The singular part of G-f-j. arises from 
the modes in the integrand with large negative energy. To 
obtain the regular part in the limit a;' — > a; we will first 
subtract this large energy contribution. For this purpose 
we define 



Vcix) 



E 

\Ei\<E^ 



fiE,)uiix)v*{x), 



(39) 



where denotes an energy cutoff introduced to regu- 
late the logarithmic divergence in i>c{x). The part of 
Vcix) dominated by the modes with large negative en- 
ergies can be approximated as [H, [H, 44] vhe{x) — 
—A{x)K{x,x;Ec) with 



K{x,x';E,) 



E 

E^<ei<E, 



2e,: 



(40) 



where here and below 4'i{x) and denote the eigenvec- 
tors and eigenvalues of the Hartree-Fock Hamiltonian 



HuF = H{n = 0)- 11 + gn{x). 



(41) 



Here fi — {^j,^ + ^4,)/2 and n{x) — {n^{x) + ni{x))/2. 
The low-energy cut-off Es in Eq. (|40|) is arbitrary, except 
that it should be chosen positive in order to avoid singu- 
larities arising from the Fermi surface. As an alternative 
to introducing a low-energy cutoff, one can add a small 
imaginary part to as done in Refs. [H, [H, In 
that case the integrand of K{x,x' ; Ec) has a peak near 
the Fermi surface, which makes the numerical integration 
over pz difficult. One can reduce this peak by increasing 
the magnitude of the imaginary part. However, that will 
worsen the large negative energy approximation of vdx). 
The low-energy cutoff which we apply here does not suf- 
fer from these problems. 

Let us next define Vsix) = i^c{x) — vhe{x). Because 
i^he{x) contains the logarithmic divergent part of vdx)^ 
the difference Us {x) is finite, and hence converges for large 
enough E^. There is some freedom in choosing vi{e{x). 



For example, we could have left out the qn(x) term in the 
Hartree-Fock Hamiltonian as in Refs. [3a,|4j]- However, 
by including this term, vhe{x) approximates Vcix) much 
better, which implies that a much smaller value of Ec is 
sufficient to compute Vs{x) accurately ji^. 

Summarizing the discussion above, we found that in 
the limit x' x, we can write Gi-ijx, t;x' , r) = Vsix) — 
A{x)K{x,x';oo). Following Refs. [23, Q the singular 
part can now be obtained by making use of the Thomas- 
Fermi approximation. In the limit x' ^ x one finds that 

K{x, x'; oo) = K{x, x; Ec') 
1 r d^k 1 



d^fc 



(2^) 
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2M 



+ iMw^p^ — ^ -f gn{p) + 



(42) 



Here 7 is an infinitesimal small positive number and Ec' 
is a second energy cutoff that can be chosen different from 
Ec- For large enough Ec' the sum of first two terms in 
Eq. (|42l) is convergent. The inhomogeneous wavevector 
cut-off kc' (p) can be found from 



+ -A/a;2p2 _ „ + gn(p) 

2M 2 ^ H'-ry yn 



Ec 



(43) 



The last term of Eq. contains the singularity. One 
can now perform the integration over k analytically, 
which in the limit x x' gives (43| , |4^ 



G^i{x,t;x',t) = - 



MA{x) 



1 



(44) 



where 



^tii^^ ''") 
A{x)M 



+ 



i^s{x) - A{x)K{x, x; Ec') 

1 . . . . ^ kc'{p) + kF{p)\ 



(45) 



Here we have introduced the length of the Fermi wavevec- 
tor kp^p) which can be found from the equation 

= M - gnip) - -Muj'p^ - i7. (46) 

To obtain A(a;) we have to multiply Eq. (P5|) by g as 
follows from Eq. (US}. 

Inserting Eq. ([B8| into Eqs. ([301) and ^ gives the 
Helniholtz free energy 



-I3\E,\ 



- J d^xGu{x,T;x,T)*A{x) 

d^xn^{x)ni{x) + fi^N^ + mNi. 



E' 



(47) 
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Although some of the individual terms in the last equa- 
tion are ultraviolet divergent, their sum, and hence the 
Helmholtz free energy, is ultraviolet finite. The diver- 
gence present in the sum over the eigenvalues of the 
Bogoliubov-de Gennes matrix is canceled by the sum over 
the eigenvalues of the Hartree-Fock Hamiltonian and by 
the logarithmic divergence originating from G^i- This 
can be made clearer by expressing G-\-i in terms of 
and K. One finds 



\Ei\<E^ 



2 



l-he 



-fi\Ei 



E 

e,<-Ec 



(48) 

In the absence of a trapping potential Ei is known ana- 
lytically; then it can be seen that Eq. (|48l) is ultraviolet 
finite. As we will see in the next section, F is also finite 
in the general case. 



V. SOLVING THE BOGOLIUBOV-DE 
GENNES-EQUATION 

In the previous section we have reduced the Dyson- 
Schwinger equation to a nonlinear equation of the form 

(A(a;),nt,^(a;),AAt,J -F(A(a;),nt,i(a^),Ait,i)- (49) 

We will now discuss how to compute the function F in 
practice. In the next section we will solve the Dyson- 
Schwinger equation. 

First we will use the symmetries of our problem to 
simplify the analysis. For zero rotation frequency the 
superfluid is in a vortex-free phase. Then, the pairing 
field A (a;) will be a function of the radial coordinate p 
only, i.e. A(a;) = A(p), where A(p) e M. 

We will assume that the first vortex that appears when 
increasing the rotation frequency carries one unit of an- 
gular momentum and is located at the center of the trap. 
The pairing field for such vortex has the following form: 
A{x) = A(p)exp(i^). After the single vortex has ap- 
peared, a vortex lattice can be formed by further increas- 
ing the rotation frequency. 

For these reasons we will make the following ansatz for 
the pairing field 



A{x) = A(p)cxp(ifc(/)), 



(50) 



where k is the winding number (unit of angular momen- 
tum) of the vortex at the center. Hence the fc = 
case corresponds to the vortex-free phase. To determine 
the onset of the vortex phase, we have to compare the 
Helmholtz free energies with A; = and k = 1. Because 
of the cylindrical symmetry of the trap and the fact that 



a possible vortex is located at the origin, the number 
densities are a function of p only, i.e. n^^^{x) — n^^i{p). 

In this case the solutions of the Bogoliubov-de Gennes 
equation, Eq. p5p . are of the following form 



Ui{x) 
Vi{x) 



1 
1 

7=L 



z _ 



1 



'2tt 



\p 



Xp), (51) 



ip.. ^^^^J^)^ (52) 



\p 



which can be verified by inserting these expressions into 
Eq. ([35)) . This also yields the Bogoliubov-de Gennes 
equation for UnmpAp) and VnmpAp)-, 



nj I Unmp^ {p) 



p I ^nmp^ (p) 

'-'nmpz \P ) 



(53) 



where 



H^{n)-p^+gnAp) A{p) 

A(p) -Hk-„i{^) + Pi-gnAp) 



and 



( d^ 1 



p2 + A2772 



2 ^ 2M 



(54) 



(55) 



Due to the factor l/y/p in Eqs. ([ST]) and (|52l) Hm{^) is 
a bit different from the single particle Hamiltonian de- 
fined in Eq. ([5]). We have inserted this factor for later 
convenience. Furthermore, by inserting the term }?rj^ 
we have modified the centrifugal potential in such a way 
that it becomes regular at p = 0. The original centrifu- 
gal potential is reproduced for 77 — 0. As we explain 
in Appendix |D1 this slight modification is necessary for 
the numerical computation of the wavefunctions and en- 
ergies. In our computations we have taken r] = 10~^ 
and checked that the results are completely stable if is 
varied by orders of magnitude around this value. 

From the normalization condition on Ui(x) and Vi{x) 
it follows that Unmp^ (p) and Vnmp^ (p) have to be normal- 
ized as 



dp [unmpApY + VnmpAp)^] = 1- (56) 



To solve the Bogoliubov-de Gennes equation numeri- 
cally, we have to discretrize Eq. ([55]) . One can do this by 
expanding the wavefunctions u{p) and w(p) in a certain 
basis. For practical purposes, this basis has to be trun- 
cated, i.e., we will represent UmnpAp) and Vmnp^(p) by 
a finite number N of basis functions ii (p) . More specifi- 
cally we will write 



N 

Ap) ^^cMp) 

i=l 



ip) 



N 



dA^{p), (57) 
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here Ci and di are the expansion coefficients, which de- 
pend on the quantum numbers n, m and p^. We will 
require the basis functions to be orthonormal in the fol- 
lowing way 



(58) 



From Eq. (I56t it then follows that the expansion coeffi- 
cients have to be normalized as 



N 

1=1 



(59) 



Equation (j53p can now be transformed into an ordinary 
eigenvalue equation for a 2N x 27V matrix which reads 



H„i — pf 
A 



A 



nmp~ 



(60) 

Here H, p and A are N x N real symmetric matrices 
which are given by 



{P'a)ij 



1 

1 

A 



dpe,{p)H„,{n)e,(p), 



(61) 



dp£,{p) [pa - gn_a{p)] ij{p), (62) 
dp£d,p)A{p)ej{p)- (63) 



Similarly we define {Hiip)ij as the representation of the 
Hartree-Fock Hamiltonian defined in Eq. (HTI) . Once 
these matrices have been computed explicitly, Eq. (pO)) 
can be solved numerically using standard linear algebra 
routines. From the solutions the wavefunctions and the 
expressions for I'six) and n^^^(a;) can be constructed. We 
will give the explicit expressions at the end of this section. 

Since truncating the basis is an approximation, we 
have to choose the basis carefully in order to make sure 
that the basis functions can describe the exact solution 
to good accuracy. We can always improve the accuracy 
of the truncation by taking a larger value of N, but the 
drawback is that this increases the computational cost as 
well. 

A good basis has to be able to describe the wave- 
functions in the case A = 0, which are the radial single 
particle wave functions given in Eq. ^ times y^. For 
that reason we will choose a basis of the following form 



e,ip) = yMp/xjk{p/x), 



(64) 



where w{x) = xexp(— a;^) and li{x) is a set of linear 
independent polynomials of degree A^ — 1. In this basis all 
single particle wave functions can be represented exactly 
by a finite number of basis functions. For the calculation 
of vortices this is a desirable feature, because in that case 
mixing between states with different angular momentum 
quantum numbers occurs. 



From Eq. ([58)) it follows that the polynomials U (x) have 
to be chosen orthogonal with respect to weight function 
w{x), i.e. 



dx w{x)li{x)lj{x) = Si. 



(65) 



The set of polynomials of increasing degree that are 
orthonormal to each other with the weight function 
x^ ex-p{—x^) on the interval [0, oo) are called Maxwell 
polynomials [i^. We will write the Maxwell polynomi- 
als for p = 1 as 4)i{x), where i = 0, 1, 2, . . . denotes the 
degree. 

One could have chosen as basis functions the Maxwell 
polynomials directly, for example li{x) ~ (j)i-i{x). How- 
ever, the computation of the Bogoliubov-de Gennes ma- 
trix, Eq. (|60p , becomes much easier if we use Lagrange in- 
terpolating functions as basis functions. These Lagrange 
interpolating functions are a particular linear combina- 
tion of Maxwell polynomials and will be specified below. 
This approach is generally known as the Discrete Vari- 
able Representation (DVR) method, or alternatively as 
the Lagrange mesh discretization [46l - l48j . By applying 
this method one can obtain hi ghly accurate values for 
the energies and wavefunctions [4^|. The DVR method 
can be applied to any set of orthogonal polynomials [HO] , 
and is not restricted to Maxwell polynomials. To our 
knowledge, this is the first time that the DVR method is 
used with Maxwell polynomials. 

The DVR method is based on the Gaussian quadrature 
formula 



dx w{x)f(x) 



N 



(66) 



Here and in the following the nodes the roots of 

the Maxwell polynomial 4in{x) and Wn the correspond- 
ing quadrature weights. The integration formula Eq. (|66|) 
is exact for all polynomials f{x) of degree less than 2N . 
From the properties of the orthogonal polynomials, one 
can show that all N roots are real and all weights are 
positive. The nodes and the weights are the only non- 
trivial properties of the Maxwell polynomials one needs 
to know in order to apply the DVR method. Since the 
Maxwell polynomials are non-standard polynomials, we 
review the computation of its nodes and weights in Ap- 
pendix [Cj 

In the DVR method one chooses the functions li{x) 
to be the Lagrange interpolating functions through the 
nodes x„. Explicitly, these functions read 



li{x) 



1 



N 

X Xji 
-, , .Xi Xr, 



(67) 



It follows directly that the polynomials li{x) satisfy the 
following useful property, li{xj) = Sij/y/wi. Because the 
combination li{x)lj{x) is a polynomial of degree 2A^ — 2, 
the Gaussian quadrature is exact for this combination. 
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Hence we can use the Gaussian quadrature to show that 
li (x) satisfies the required orthonormality condition given 
in Eq. ([eS]). 



dx w{x)li{x)lj (x) 



N 



Wk 



(68) 



We will now compute the matrices that appear in the 
Bogoliubov-de Gennes equation, Eqs. (|6T1) - (|63|) explicitly. 
In the DVR method one uses the Gaussian quadrature 
to compute the integrals. In this way we find that 

{p-a)ij ~ ^aSij - gn^a{Xxi)dij, (69) 
(A)y « A(Ax05^,. (70) 

Here one of the nice features of the DVR method appears. 
All parts of the Bogoliubov-de Gennes matrix that con- 
tain no derivatives become diagonal and are easy and fast 
to evaluate. Furthermore, the values of A(p) and n-f-4(p) 
have to be evaluated only at the mesh points p = Xxi. 
These mesh points are unevenly spaced. 

Another attractive feature of the DVR method is that 
matrix elements of derivative operators can be com- 
puted exactly. We discuss the computation of H in 
Appendix |D] The exact expression for the matrix H 
can be obtained by inserting Eqs. (ID14I) and (|D17|) into 
Eq. (jDip . While H becomes much more complicated 
than fi and A, this is not a disadvantage, because we 
only need to compute H once. This is in contrast to /2 
and A which will change from iteration to iteration when 
solving the Dyson-Schwinger equations. 

In order to describe a single particle wave function with 
quantum numbers n and m in this basis exactly, we need 
to take N > 2n + \m\. For a given value of N therefore 
only the lowest {N — \m\)/2 eigenvalues and correspond- 
ing eigenvectors can be computed exactly in this basis. 

Let us now order the eigenvalues of the Bogoliubov-de 
Gennes matrix and the Hartree-Fock Hamiltonian, i.e. 
label them in such a way that 



■«+i,mp,, n = l. ..2N, (71) 
< en+i,nip,, n = l...N. (72) 



If we can compute the first n,„ax eigenvalues of the sin- 
gle particle Hamiltonian accurately then the eigenvalues 
of the Bogoliubov-de Gennes matrix with Ui < n < n{ 



where rii — N 



1 and rif — A^ + rimax, will be com- 



puted accurately. We will write the maximum n quantum 
number as 



kN — max(|m|, \k — m|) 



(73) 



where [x\ denotes the fioor function, and k is a free pa- 
rameter which equals 1 if all single particle eigenvalues 
and vectors need to be computed exactly. Typically the 
next few eigenvalues and eigenvectors can also be com- 
puted very reliably, although they deteriorate rapidly at 
some point. A larger value of k is advantageous because 



more eigenvalues and vectors are taken into account in 
the same basis. From our experience one can safely take 
1 < K < 1.2. The angular quantum number is varying 
between mmin and mmax- These are in the case of fc > 
given by 



-[nN\+k, 
IkN\. 



(74) 
(75) 



To solve the Bogoliubov-de Gennes equation we only 
need to evaluate Vs(,p) and n^,],ip) at p = Xxi. At the 
mesh points the wave functions become rather simple and 
read 



^nmpz {,XXi ) 
^nmpz {Xxi ) 




exp(-x2/2), (76) 
exp(-x2/2). (77) 



Combining now the last two equations with Eqs. (1511) and 
([5^ we find that at the mesh points the densities read 



n^{Xxi) 
niiXxi) 



27rA^w, 



J E fi^nmpjcl (78) 



n— ni 



2nX^w, 

Here we introduced the symbol 



-EnmpMl (79) 



(80) 



where pc is a cutoff on the pz integration. Furthermore, 
we used the fact that all integrands are symmetric in pz- 
We have performed the integration over pz numerically 
using the adaptive Simpson method. After that we per- 
formed the sum over m. 

The function I'six) from which A(a;) can be obtained 
becomes at the mesh points 



+ A(AxO ^ - 



Es) 



(81) 



where €nmp^ and Ci are respectively the eigenvalues and 
eigenvectors of the matrix representation of the Hartree- 
Fock Hamiltonian, {HY^p)ij- To compute A(Axi) we 
also need to evaluate K{Xxi,Xxi,Ec'), which is defined 
in Eq. pO)) . To obtain K only the eigenvalues and 
eigenvectors of the Hartree-Fock Hamiltonian in the case 
Pz = have to be computed numerically. Their values for 



nonzero pz follow trivially, since 



,o+P?/(2M). 
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We then obtain 



27rA2 



Wi 



j;2 mmax ^Imax 

m— mmin n— 1 

dp2_ 1 



2TTh 2e 



nmpz 



(82) 



The value of Ec' is hmited by the largest enmp^ one 
can compute reliably. An estimate of this value is 
Ec' < kN — fi + gn{0). The integration over in 
Eq. (15^ can straightforwardly be performed analytically 
(we will however not write down the result here). There- 
fore K{Xxi, Xxi] Ec') can be computed much faster than 
VsiXxi) for which numerical integration over is re- 
quired. To compute K{Xxi, Xxf, Ec') we can thus afford 
to use a much larger basis (we have taken N = 128) than 
for i^six). A larger basis implies that we can choose a 
larger value of the cutoff Ec' leading to more reliable an- 
swers. The values of -fC on a coarser mesh can then be 
obtained from interpolation. 

From Eq. (pS)) we obtain the pairing field 

A(Axj) = gVsiXx^) ~ gA{Xxi)K{Xxi, Xxi] Ec') 
A{Xx,)M 



+9- 



1 



kc'iXx,) (83) 
kc'{Xxi) + kpiXxi) ' 



2''Fi^^^)^o^ykc'{Xx,)-kp{Xx,) 



The Helmholtz free energy per unit of harmonic oscil- 
lator length in the z-direction follows from Eq. p8)) and 
reads 



T = -X 




2TrX^ y^^WiVsjXxi)* A{Xxi)e'^^i 



i=l 
N 



g2TiX^ ^ w^n^{Xxi)n^{Xxi)e^- + fi^Af-^ + mJ^i- 

(84) 



Here we have used the Gauss-Maxwell quadrature to 
compute the integrals over p. It can be shown numer- 
ically that the integrand of T decreases rapidly for large 
Pz, making T ultraviolet finite. 

Another important quantity which now can be con- 
structed is the expectation value of the angular momen- 
tum operator L^- With help of the number densities, 
Eqs. (1751) and ([75]) . the angular momentum density can 



be written as 
/ 



n— rii 

+ ik-m)dff{-EnmpJ]- (85) 



Integrating the last equation over p and (j) gives the an- 
gular momentum per unit harmonic oscillator length in 
the z-direction which we will denote by C^. 



VI. SOLVING THE DYSON-SCHWINGER 
EQUATION 

Solving the Dyson-Schwinger equation amounts to 
finding the solution of Eqs. ([751), dZll) and dHS]), together 
with the constraint on the number of particles. Schemat- 
ically the equation to be solved is of the form given in 
Eq. Such equation can be solved using a multi- 

dimensional root-finding method. For that purpose we 
will use the Newton-Broyden method [H^, which leads 
to very fast convergence once close to the solution. As 
an input to the Newton-Broyden method one should pro- 
vide initial guesses for A, n-f-^, and and also provide 
the Jacobian of F. 

To obtain the initial conditions we will use the 
Thomas-Fermi approximation. The Thomas-Fermi ap- 
proximation to the density is found by solving the fol- 
lowing equation 

(p) ^ / ^Pf(^ + -Muj^p" - QLz 
J {27rh)^-' ^2M 2 ' 

- Mt4 +.9%,t(p))- (86) 

If r = the last equation becomes 



n-\,i{p) 



1 



67r2A3 



a;2 ; A2 



3/2 



(87) 



This equation can easily be solved numerically. If /i-f- = 
the Thomas-Fermi radius (the minimal p for which 
n^_^(p) — 0) of the gas becomes 



R = 



1 -fi2/tj2 




The initial values for can be obtained by integrating 
the Thomas-Fermi density profiles and solving numeri- 
cally for the desired A/f^j,. 

For the initial condition to the pairing field in the case 
fc = we use the result of the BCS theory in the weak 
coupling limit (see e.g. Ref. (Hj). 



A(p) =4fcF(p) A-^cxp -2- 



2kF{p)\a\ 



fujj. 



(89) 



where kp{p) is defined in Eq. (j46p . As an initial condition 
for nonzero fc, we multiply this equation by a factor 1 — 
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exp[— where the BCS coherence length equals 



fcF(0)A^ 
7rA(0)/(fix^)' 



(90) 



In the Newton-Broyden method the Jacobian has to 
be computed at each iteration. We have computed the 
initial Jacobian using finite difi^erences. If /i-f- = /ij, there 
are 2{2N + 1) evaluations of F required, so obtaining the 
initial Jacobian is computationally expensive. 

However, it is not necessary to apply finite differences 
to recompute the Jacobian in the next iteration. The 
next Jacobian can be obtained from the previous one us- 
ing a Broyden update [56|. Such update has negligible 
computational cost. After the initial Jacobian has been 
acquired Newton-Broyden iterations typically reach con- 
vergence in about ten iterations. The time it takes to 
reach convergence is completely dominated by the time 
it takes to compute the initial Jacobian. The computa- 
tional cost of the whole problem scales roughly with N^, 
a factor of originates from solving the Bogoliubov-de 
Gennes equation, a factor of about N arising from the 
sum over angular quantum numbers and another factor 
of N from the Jacobian. 

Given the solution, the Helmholtz free energy can 
straightforwardly computed by applying Eq. (|84l) . 



VII. RESULTS 

In the numerical computations the following param- 
eters were used rj = 10^^, k = 1.1, Eg — 2A3hu!, 
Pc/fi — 20/A and Ec' = lOOhoj. The adaptive Simp- 
son integration was performed with a relative precision 
goal of 10~^ and absolute precision goal 10~^° for each 
individual component of A{Xxi) and n-f,i{Xxi). The ac- 
curacy goals for the computation of the initial Jacobian 
were set much lower, which speeds up the computation. 
As long as these goals are not set too low they will not 
ruin the convergence of the Newton-Broyden algorithm. 
We would like to stress that a less accurate Jacobian that 
can make the Newton-Broyden algorithm converge does 
not infiuence the accuracy of the final values of A{Xxi) 
and n^^i{Xxi). 

The iterations of the Newton-Broyden algorithm 
were performed until the relative difference be- 
tween the norm of {A(Xxi),n-f^i{Xxi),Aff^i) and 
F{A{Xxi),n-f^l{Xxi),fi^_l) became less than 10^''. 
Frequently, a relative accuracy of 10~^ could be reached. 

The number of basis functions N was chosen such that 
convergence was reached. At least all energy levels below 
the Fermi energy have to be computed accurately. Since a 
larger number of particles implies a larger Fermi energy, 
for a larger number of particles a larger value of N is 
required. We have done computations from = 16 to 
A^ = 80. A calculation with A^ = 80 took several days 
on a single modern CPU. For calculations with Af = 100, 
200 and 1000 we typically used respectively A^ = 40, 48 
and 64. 



fcF(0)|a| 




FIG. 3: Inverse interaction strength at the center of the trap 
as a function of scattering length, for f2 = and different 
number densities of particles A/". 



We have checked that our results are completely sta- 
ble under acceptable variation of these parameters. By 
studying these variations, we have convinced ourselves 
that the largest values of A{Xxi), n-^.i{Xxi) have a rela- 
tive accuracy of a least 10"'^. The free energy could be 
obtained with a relative accuracy of 10~^. 

Furthermore we have taken T = throughout and con- 
sidered the situation that the number of particles per 
unit harmonic oscillator length in the z direction in each 
species is equal, i.e. J\f-f — A/'j. = J\f/2. We have investi- 
gated situations with J\f = 100, 200 and 1000, different 
scattering lengths and different rotation frequencies. The 
scattering lengths we have considered correspond to in- 
verse interaction strengths at the center of the trap in the 
range 1 < l/(fcF(0)|a|) < 4. We have depicted their rela- 
tionship at zero rotation frequency in Fig. [31 In this range 
the Hartree-Fock-Bogoliubov approximation is expected 
to be valid. If one wants to study stronger interactions 
one has to take into account the higher order diagrams 
in order to get a reliable result. 

To get an idea of the scales in a typical experiment, 
we can use that in Ref. [l| a Fermi gas made out of ^Li 
atoms was studied in a trapping potential with radial 
frequency w/(27r) = 57 Hz. This situation corresponds 
to A ^ 5.4 /im and huj/ks ^ 2.7 nK. 

We will now first give a detailed overview of the results 
at zero rotation frequency. After that we will discuss the 
effects of rotating the trap and present the main objective 
of this work, the critical rotation frequency for vortex 
formation. 



A. Zero rotation frequency 

In Fig. m we compare the pairing field of the vortex- free 
superfiuid (k = 0) with the pairing field of the vortex 
with unit angular momentum (k = 1). These pairing 
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FIG. 4: Pairing field as a function of radius, for a — — O.IA, 
Af = 1000 and = 0. 




FIG. 5: Number density as a function of radius, for a — 
-O.IA, AA= 1000 and n = 0. 

fields were computed for A/" — 1000 and a — —O.IA with 
N = 64. Because the data points lie so close to each 
other, we have only displayed a line that interpolates 
through the data points for visibility reasons. Since the 
rotation frequency was taken to be zero, the vortex is 
metastable. 

The pairing field clearly vanishes at the center of the 
vortex. Away from the vortex core the pairing field is 
restored to its value in the fc = case. The typical dis- 
tance at which this happens is the BCS coherence length 
^. This implies that the size of the vortex core grows 
when decreasing the strength of the interaction. For very 
weak interactions, ^ can become larger than the radius 
of the gas. In that case even a metastable vortex is no 
longer possible. In the second part of this section we will 
study the size of the vortex core at the critical rotation 
frequency for vortex formation in some detail. 

The vortex also leaves its imprint on the corresponding 
density profiles which are displayed in Fig. [Sj As already 
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FIG. 6: Pairing field as a function of radius, for M — 100, 
k — and different scattering lengths. 

found in Refs. [ill, [13]; the density at the center of the 
trap is significantly depleted in the presence of a vortex. 
To compensate for the removal of particles at the center, 
the gas will expand. This is a tiny effect and is, therefore, 
not visible in Fig. [5l In the next subsection we will study 
density depletion at the critical rotation frequency. 

The density profile for normal pairing is very well de- 
scribed by the Thomas-Fermi approximation, which is 
the solution of Eq. (|57|) . The pairing field however only 
agrees qualitatively with the Tomas-Fermi approxima- 
tion, Eq. ([55]) . as was also concluded in Ref. (36j . 

We have displayed pairing field profiles for Af — 100 in 
Fig. ini The interaction strength was taken to be weak, in 
the range 0.07 < |a|/A < 0.09, which leads to small pair- 
ing fields. In such situations we encountered oscillations 
in the pairing field. To ensure that this is not a numerical 
artifact, we have compared these profiles computed with 
AT = 32, iV = 48 and N = 64. We find that they are 
completely consistent with each other. The oscillations 
in the pairing field are only prominent in the case of a 
small number of particles with weak interactions. 

Let us now study the effect of variation of the scatter- 
ing length and the number of particles in the zero rota- 
tion limit. In Fig. [7] we display the pairing field at the 
center of the trap as a function of the scattering length. 
Clearly, increasing the interaction strength and the num- 
ber of particles (i.e. the Fermi wave number) both lead 
to larger pairing fields. This is qualitatively in agreement 
with the BCS pairing formula, Eq. ([55]) . 

In Fig. [5] we have displayed the corresponding number 
density at the center of the trap. Increasing the number 
of particles leads naturally to a larger number density 
at the center. Stronger attractive interactions lead to a 
more compressed gas, which likewise results in a larger 
density at the center. 

In Fig. [S] we have displayed the corresponding chemi- 
cal potential. Obviously a larger number of particles im- 
plies a larger chemical potential. The chemical potential 
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FIG. 7: Pairing field at the center of tlie trap, as function 
of the scattering length, for different number of particles and 
n = 0. 
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FIG. 8: Density at the center of the trap, as function of the 
scattering length, for different number of particles, and = 0. 
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FIG. 9: Chemical potential as function of the scattering 
length, for different number of particles, and f2 = 0. 



FIG. 10: Pairing field, number density, and angular momen- 
tum density as a function of radius. The results correspond 
to a vortex-free superfluid with A/" — 1000, a — — O.IA, and 
Q, = O.lSoj. 
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FIG. 11: Same as in Fig. 1101 but for a vortex with k = 1. 



decreases with increasing the strength of the interaction. 
This is because the Hartree term —gni^-f(p), which acts as 
a sort of inhomogeneous chemical potential, grows with 
increasing the interaction strength. 

The radius of the gas can, to very good approxima- 
tion, be obtained by inserting the values of the chemical 
potential in the Thomas-Fermi estimate, Eq. (|88)) . This 



gives at zero rotation frequency R/X = y^2fi/huj. It then 
follows from Fig. ^ that increasing the number of parti- 
cles increases the radius. On the other hand, increasing 
the interaction strength reduces the radius. The gas be- 
comes more compressed because the interaction between 
the two components is attractive. 
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B. Non-zero rotation frequency 

To illustrate the effects of rotation, we have displayed 
the pairing field, the number density and the angular 
momentum density for f2 = O.lSw, Af = 1000, and a = 
-0.1 A in Fig. [TO] (A: = 0) and Fig. [lU (fc = 1). These 
figures can be compared to the results at zero rotation 
frequency which are displayed in Figs. |4] and [H As we will 
show below, = 0.15w is the critical rotation frequency 
for vortex formation in this situation. 

A vortex-free superfluid state cannot carry angular mo- 
mentum. For that reason the angular momentum density 
vanishes in the region where the pairing field is sizable, as 
can be seen in Fig. 1101 Since angular momentum density 
appears at large p, it indicates the presence of unpaired 
fermions at the edges of the gas. Another way to ob- 
serve this effect is that for large radial coordinates, the 
pairing field disappears before the number density does. 
It can be seen in Fig. [Tl] that a vortex generates angu- 
lar momentum density in the superfluid region. For the 
same reasons as in the fc = case, unpaired fermions are 
present at the boundaries of the gas. 

By careful comparison of Figs. fTOl and [Til one can ob- 
serve that the pairing field of the vortex is slightly larger 
in the outer region. This is generally the case and leads 
in addition to the effects mentioned in the introduction 
to a fourth contribution to the energy difference between 
a vortex and a vortex-free phase. Rotation increases the 
radius of the cloud as well. However, at this rotation rate 
this is only a very small effect and is therefore not visible 
in the figures. 

Now let us discuss the determination of the critical ro- 
tation frequencies for unpairing and vortex formation. To 
obtain these frequencies we have computed the Helmholtz 
free energy. The phase with the lowest free energy is the 
preferred phase. 

In Fig. [TO] we have displayed the Helmholtz free en- 
ergy divided by the number of particles, for A/" = 1000 
and a = —0.1 A. A number of interesting features of the 
superfluid are shown in this figure. First of all, the super- 
fluid phase is always preferred over the unpaired phase, 
since A = has the largest free energy. Furthermore, 
for < 0.05ci; the gas forms a vortex-free superfluid. 
It can be seen that in this region the free energy does 
not depend on the rotation frequency. This indicates 
that the entire gas is in a superfluid state. However, for 
il > 0.05a; the free energy of the vortex-free phase starts 
to decrease when increasing the rotation frequency. This 
implies that the gas has acquired angular momentum, 
which occurs via unpairing the fermions at the edges of 
the gas. Hence for O.OSw < il < 0.15w the gas forms 
a vortex-free superfluid with unpaired fermions at the 
boundaries. At il > 0.15a; a superfluid with a, k — 1 vor- 
tex becomes the preferred phase. One can see this more 
clearly in Fig. [TO] where we have displayed the difference 
in free energy between the vortex phase with fc = 1 and 
the vortex-free phase. The critical rotation frequency can 
be found from interpolation of the data points which in 
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FIG. 12: Helmholtz free energy divided by the number of 
particles, as a function of rotation frequency, for a = — O.IA 
and M — 1000. The label fc = corresponds to a superfluid 
without vortices, the nonzero values of fc correspond to a sin- 
gle vortex at the center of the trap with angular momentum 
fc. The label A = corresponds to the situation in which all 
fermions are unpaired. 
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FIG. 13: Difference in Helmholtz free energy per unit length 
in the 2-direction between the vortex phase with fc = 1 and 
the vortex-free phase with fc = 0, as a function of rotation 
frequency, for a — O.IA and TV — 1000. 



this case yields ilc — 0.149a;. For il < 0.15a; the fc = 1 
phase is metastable. Because superfluids with a vortex 
carry angular momentum, the derivative of their free en- 
ergy with respect to rotation frequency is negative, even 
at zero rotation frequency. In the unpaired phase this 
derivative vanishes at zero frequency, because the fully 
unpaired gas does not contain angular momentum at zero 
frequency. 

At zero temperature, one can also compute the critical 
rotation frequency for unpairing in a more direct way. At 
zero rotation frequency all quasi-particle excitations (ex- 
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FIG. 14: Phase diagram of a two-component Fermi gas as 
a function of scattering length and rotation frequency, for 
M — 1000 and T = 0. The hues correspond to the phase 
boundaries. The label A indicates that the entire gas is in a 
vortex-free superfiuid state. The label B indicates a vortex- 
free superfluid in the center with unpaired fermions in the 
outer regions of the gas. The label C indicates a superfluid 
with vortices in the center and unpaired fermions in the outer 
regions. 



cept the superfluid phonon) are gapped, i.e. \Enmp, \ > 0. 
As follows from the discussion in Appendix [B] if /^t = /^i 
and = 0, both Enmp^ and —Enmp^ are eigenvalues of 
the Bogolibuov-de Gennes matrix. Rotation shifts these 
eigenvalues downwards by mMl. As long as no gapless 
mode arises the rotational contributions from positive 
and negative energies cancel so that this shift has no 
effect on the free energy. For that reason the free energy 
for k — Q stays constant up to a certain rotation fre- 
quency. Only when the first gapless mode appears, the 
free energy will change. The minimal rotation frequency 
at which this occurs is the critical rotation frequency for 
unpairing, fi^. Thus this rotation frequency can be found 
from the solutions at = in the following way 



1 



(91) 



where the minimum is to be taken over all values of n, 
TO, and Pz- Determination of VL^ in this way is computa- 
tionally much less expensive than obtaining it from the 
free energy. 

Let us now discuss our main result: the phase diagram 
as a function of scattering length and rotation frequency. 
In Figs. [T4l and [T5] we have displayed these diagrams for 
TV = 1000 and TV = 200 respectively. 

There are two transitions in these phase diagrams. The 
lower line denotes the unpairing transition. The order 
parameter corresponding to this transition is the angu- 
lar momentum. This transition is of second order since 
the angular momentum changes continuously. At T 7^ 
this transition turns into a crossover. Hence at f2 = fit. 
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FIG. 15: Same as in Fig. [H but for M = 200. 



and T — {) the gas resides at a so-called quantum criti- 
cal point. Above this critical point the order parameter 
behaves as Cz ~ where t = — 0„)/J7„. We find nu- 
merically that the critical exponent has the value (3 = 1. 

The upper line denotes the critical rotation frequency 
for the formation of a vortex with unit angular momen- 
tum at the center of the trap. Energy arguments suggest 
that with increasing rotation frequency the first vortex 
configuration that will nucleate is a single vortex with 
fc = 1. A single vortex with fc > 1 will have lager energy, 
as can be seen in Fig. [T^l Several vortices with fc = 1 
have again larger energy and their nucleation would re- 
quire larger than critical rotation frequencies. Therefore, 
the upper line shows the critical rotation frequency for 
vortex formation. The order parameter corresponding 
to this transition is the winding number of the vortex. 
Since this winding number changes discontinuously, this 
transition is of first order. 

Increasing the absolute value of the scattering length 
leads to a larger critical rotation frequency for unpairing. 
Furthermore, for a given scattering length f2„ becomes 
larger when the number of particles is increased. Both 
effects can be explained by the fact that a stronger bound 
pair is more difficult to break. 

As can be seen from the phase diagrams, we find that 
vortices are formed only for relatively large negative scat- 
tering lengths. The critical rotation frequency for vortex 
formation has a minimum at a certain intermediate value 
of the scattering length. This minimum arises from the 
interplay of two effects. Firstly, the energy cost of creat- 
ing a vortex at zero rotation frequency increases with in- 
creasing the negative scattering length. This explains the 
rise of the critical frequency at large negative scattering 
lengths. The second effect is caused by the difference in 
energy gain due to rotation. The fc = 1 phase will always 
have a larger rotational energy gain than the fc = phase 
due to the angular momentum generated by the vortex. 
However, for small interaction strengths above the un- 
pairing transition, the difference between these gains is 
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FIG. 16: Relative density depletion at the core of a vortex, 
as a function of scattering length, at the critical rotation fre- 
quency for vortex formation. 
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FIG. 17: Half-width of the vortex as a function of scattering 
length, at the critical rotation frequency for vortex formation. 



relatively small. This is because in this case it is rela- 
tively easy to break the pairs at the boundaries of the 
gas, which contribute to the rotational energy gain in 
both the k = and fc = 1 phase. As a result of this 
effect the critical frequency increases for small negative 
scattering lengths. For a certain small scattering length 
the difference in rotational energy gain cannot overcome 
the costs associated to the vortex. For this reason at 
small negative scattering lengths the vortex phase has 
an abrupt transition to a vortex-free phase. We see that 
below a « -0.155A for M = 200 and a « -0.085A for 
Af = 1000 vortex formation does not occur at all the 
rotation frequencies displayed in the phase diagram. 

Vortex formation sets in at a lower rotation rate when 
the number density of particles is increased from Af = 200 
to J\f = 1000. For the number of particles we have inves- 
tigated we find that the vortices always appear together 
with unpaired fermions at the edges of the gas. One 
could speculate that for a larger number of particles ftc 
will be reduced so that a vortex phase will appear before 
unpairing at the edges could become possible. In other 
words, we anticipate that for a large number of particles 
the phase diagram might feature a direct phase transi- 
tion from the A to the C phase without the intermediate 
B phase. 

In Fig. we have displayed the central number den- 
sity of a vortex over the central density of a vortex-free 
superfluid, at the transition to vortex formation. It can 
be seen that the amount of density depletion is relatively 
small for weak interactions and grows with increasing the 
interaction strength. 

Let us define the half-width d of the vortex to be the 
radius at which A{p)k=i/ A{p)k=o = 1/2. In Fig. [l7] 
we have displayed this half-width at the transition for 
vortex-formation. For M = 1000 it can be clearly seen 
that weak interactions lead to larger vortices, which is 
caused by the increase of the BCS coherence length. 



VIII. CONCLUSIONS 

In this article we studied a two-component Fermi gas 
with attractive s-wave interactions confined in a cylin- 
drically symmetric harmonic trap. Our key results are 
summarized in a phase diagram spanned by the rotation 
frequency and scattering length for zero temperature and 
a fixed number of particles. Explicit results are shown for 
a number density of 1000 and 200 particles per unit har- 
monic oscillator length in the z-direction in Figs. [HI and 
[T5] respectively. 

To obtain the phase diagram we have used the two- 
particle irreducible effective action. We only took into ac- 
count the leading order diagrams, which is equivalent to 
the Hartree-Fock-Bogoliubov approximation. This con- 
strains our study to interaction strengths of magnitude 
l/(fcj?(0)|a|) < 1. The equations we obtained were solved 
numerically using the DVR method based on Maxwell 
polynomials. 

In the phase diagram three phases can be distin- 
guished. For small rotation frequencies the entire gas 
forms a superfiuid. At a certain critical frequency a sec- 
ond order transition occurs to a superfluid phase, which 
features unpaired fermions that are concentrated at the 
edges of the gas. At this critical rotation frequency the 
gas resides at a quantum critical point when the tem- 
perature vanishes. For even larger rotation frequencies 
vortices are formed via a first order transition. These 
vortices only appear for large negative scattering lengths. 
We have found that at a certain scattering length the crit- 
ical rotation frequency for vortex formation has a mini- 
mum. 

The presence of unpaired fermions at the boundaries of 
the gas results in an increase of the critical rotation fre- 
quency for vortex formation. For this reason one cannot 
use the free energy difference between a vortex phase and 
vortex-free phase at zero rotation frequency to compute 
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the critical rotation frequency for vortex formation. 

Our theoretical findings can be compared to the ex- 
periment at will, since the vortices have been observed in 
rotating two-component Fermi gases [HQ- It would be 
interesting to obtain the structure of the phase diagram 
from such experiments. 

The theoretical understanding of the phase diagram 
can be improved in several ways. It would be worth- 
while to investigate how the phase diagram is modified 
by temperature and by an imbalance in the number of 
fermions. Furthermore, it would be very useful to extend 
our analysis to larger interaction strengths, in order to 
obtain reliable results in the unitary regime. 

Finally we would like to point out that vortices can 
also be induced by synthetic magnetic fields, as has been 
shown experimentally in a Bose-Einstein condensate [5l| . 
Since such synthetic magnetic field is similar to rotation, 
another interesting extension of our work would be to 
compute the critical synthetic magnetic field strength for 
vortex formation in a two-component Fermi gas. 
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Appendix A: The 2PI effective action 

Consider an action S — — ?i5'^Gq -I- ^int, where Gq ^ 
denotes the bare inverse propagator and ^int is the part 
of the action that contains interactions. The partition 
function Z[K] corresponding to this action in the pres- 
ence of a source term K is given by 

Z[K]= J v¥V'^exp[-S/h+ . (Al) 

Let us denote the exact propagator in the presence of a 
source term K as G. The 2PI effective action is defined 
as [13 

r[G] = - log Z[K] + Tr(XG), (A2) 

where K has to be chosen in such a way that exact prop- 
agator in the presence of K equals G. Since the exact 
propagator has to satisfy the Dyson-Schwinger equation, 
we can conclude that K = G'^ -G^^ + Ts[G], where E[G] 
is the IPI self-energy. Taking the derivative of Eq. (|A2p 

with respect to G gives = K, so that in the ex- 

tremal points r[G] = — logZ[0]. Inserting the expression 
for K in Eq. ((X2|) gives 

r[G] = - logZa - Tr (Go ^G - 1) + Tr (S[G]G), (A3) 



where Z2 = Z[G-'^ - G^'^ + T,[G]], which implies that 
Z2 is the partition function of a theory with action S2 — 
-h^^G-^^ + S2Mt, where 52,int = -;i*^S[G]* + S^t- 
One can now compute — log Z2 in a pcrturbative series, 
using G as the propagator. It follows that the first two 
terms in this pcrturbative series are 

- logZz = -TrlogG^i - Tr(S][G]G) + . . . . (A4) 

Because the IPI self-energy is now included in the inter- 
action term, cancellations will occur such that only the 
2PI diagrams survive [13 ■ Hence 

r[G] = -TrlogG-i -Tr(GoiG- 1) +r2[G], (A5) 

where r2 [G] is now the sum of all 2PI diagrams generated 
by the interaction S'int with propagator G. 

Appendix B: Derivation of the Bogoliubov-de 
Gennes-equation 

The inverse of a non-singular Hcrmitian matrix A can 
be obtained from its eigenvalues A„ and the correspond- 
ing orthonormal eigenvectors |n), with {n\m) = 5nm- 
Putting the n-th eigenvector in the n-th column of a new 
matrix J7, one finds that U is unitary, i.e. U^U ~ 1 and 
A — UhU\ where A = diag(Ai, A2, . . .). The inverse 
of A can now be constructed as A^^ = C/A^^C/'^, where 
A^^ ~ diag(l/Ai, I/A2, . . .). By performing the matrix 
multiplications, the last equation can be conveniently 
written as A^^ = \n){n\/\n- A single component 
of the inverse matrix now reads A~^^ = ^n^i\'n) {n\j) / Xn- 

The inverse Nambu-Gor'kov propagator G~^, Eq. p3|) . 
can be written as —hG~^ — hd/dr + "H, where the Hcr- 
mitian matrix "H is given by Eq. (j34p . 

Since H is independent of r, the eigenfunctions of G~^ 
are a product of eigenfunctions of hd/dr and eigenfunc- 
tions of H. Because G{X, X') has to satisfy anti-periodic 
boundary conditions in imaginary time, the properly nor- 
malized eigenfunctions of hd/dr are plane waves 0m (t) = 
exp(— iti;,„r)/-y/H^ with eigenvalue — iftwm, where the 
Matsubara frequency = (2m-t- l)7r/(/i/3), with m S Z. 
Let us denote the normalized eigenfunctions of % as 
{ui{x),Vi{x))^ with corresponding eigenvalue Ei (which 
is real). This eigenvalue equation which is given explic- 
itly in Eq. (I55|) is known as the Bogoliubov-de Gennes 
equation |42|]. Normalization (t/^C/ = 1) implies that 

j <l\[\u,{x)\^ + \v,{x)\^]=l. (Bl) 

Furthermore from completeness {UW = 1) one finds 

( u,ix)u*ix') u,{x)v*ix') \ ^ 
y v^{x)u*{x') v,{x)v*{x') J 

5{x~x')(^l ly (B2) 



18 



We can now invert the inverse propagator G ^ , to ob- 
tain the Nambu-Gor'kov propagator, 



G{X,X') 



1 



■EE 



1 



Ui{x)u*{x') Ui{x)v*{x') 
v,{x)u*{x') Vi{x)v*{x') 



(B3) 



We can see from this equation exphcitly that 
G^^{X,X') = G^i{X',X)*. In order to obtain the pair- 
ing field and the number densities we need to evaluate 
G{X, X'^). Here X'^ = (cc, r ± 77) with 77 an infinitesimal 
small positive number. In this limit one can compute the 
sum over Matsubara frequencies exactly. After using the 
completeness relation, Eq. (IB2p . one then finds 



G{X,X'^) 



STffp)! Uiix)u*{x') u^{x)v*{x') 
2^^^^''\ v,{x)u*ix') v,ix)v*{x') 



-9{t)5{x-x') 



1 
1 



(B4) 



where f{E) = [exp{f3E) + 1]^^ denotes the Fermi-Dirac 
distribution function and 6{x) is the unit-step function. 
The term proportional to the step function refiects the 
anti-commutation relation of the fermionic operators. 
The sum over i runs over all eigenvalues. Using Eqs. (|23p. 
(PSj). and (|26p one can now read off the expressions for 
the pairing field and the number densities. 

If /i-]- — ^xi, then the densities of the two species are 
equal so that nf{x) = n^_{x). By taking the complex 
conjugate of Eq. (1551) it follows in this case that if Ei 
is an eigenvalue of H with eigenvector {ui{x),Vi{x)) , 
then also ^Ei is an eigenvalue of Ji with eigenvector 
{v*{x),—u*{xj)'^. One can now use this fact to restrict 
the sum over n to eigenvectors with positive eigenvalues 
only, so that if /i-f- = /ij. one has 



G{X,X'^)^ J2 fiEr 

Ei>0 

+ ^ [!-/(£;,)] 

Ei>0 



u,{x)u*{x') u,{x)v*{x') 
v,{x)u*{x') v,{x)v*{x') 



v*{x)vi(x') -v*{x)ui{x') 
-u*{x)v^{x') u*{x)ui{x') 



9{T)Six-x') 



1 
1 



(B5) 



From this equation one can read off the expressions for 
the number density and pairing field as they often appear 
in the literature. In this paper however we will solely use 
Eq. (jB4l) . because it leads to more compact expressions 
and has broader validity. The only slight disadvantage is 
that in Eq. (|B4p we have to sum over all eigenvalues. 

To evaluate the grand potential, Eq. ([HO)) , we need to 
compute TrlogG"^. Using that the trace of a logarithm 
is the sum over the logarithm of the eigenvalues one finds 

iTrlogG-i = ^ E E log(-i^^™ + (B6) 



In order to perform the sum over the Matsubara frequen- 
cies one adds and subtracts the following infinite constant 
to the last equation 



log(-iSw,„) - log(2) 



(B7) 



After summing over Matsubara frequencies one finds 



TrlogG" 



E 



\E^ 
2 



+ ^log(l + e-^l^-l) 



C. 



(B8) 

Since C is independent of Ei it shifts the thermodynamic 
potential by an irrelevant constant and can therefore be 
ignored. Now the result Eq. (|B8I) is not entirely correct. 
For example it is still infinite and in the limit of A (a;) = 
the grand potential of an unpaired Fermi gas is not ob- 
tained. To cure this problem one needs to take carefully 
the limit 77 — 0. We proceed as in Refs. [s^,!!^ to obtain 



TrlogG" 



E 



\E^ 
2 



log 1 



E< 



(B9) 



here are the eigenvalues of the Hartree-Fock Hamil- 
tonian, which is defined in Eq. (HI]). The last equation 
is indeed finite and reduces to the grand potential of an 
unpaired Fermi gas in the case A(a;) = 0. 



Appendix C: Computation of nodes and weights of 
Maxwell polynomials 

Any set of orthonormal polynomials of increasing de- 
gree i and hence also the Maxwell polynomials {x) sat- 
isfies the following recursion relation (see e.g. Refs. [H, 
0), 



\/ /3i+i4)t+i{x) ^ (x ~ ai)4}i{x) - y/l3i4)i-i{x). (CI) 

Here a,; and /3i are the recursion coefficients. 

Once the recursion coefficients are known, the nodes 
(which are the roots of 4)m{x)) and weights of order N 
can be found by solving the following eigenvalue equation 
numerically (see e.g. Refs. [53, 54]), 



\ ( '/>o(a;„) ^ 

(l>l{Xn) 



( MXn) \ 

<^i(a;„) 

(j)N -liXn) J 



(C2) 
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The eigenvalues of the last equation are the A'' nodes a;„. 
The weights follow directly from the eigenvectors in the 
following way (see e.g. Refs. [13, [131), 



w„ 



'N-l 

E 

. i=0 



(C3) 



The eigenvectors have to be normalized in such a way 
that that the orthonormality condition is satisfied, so 
that (j^oixn) — [J^ dxw{x)]^^/^. 

The recursion coefficients can be computed using the 
Stieltjes procedure, see e.g. Ref. [55|. Although the re- 
cursion coefficients for the Maxwell polynomials can be 
computed analytically in this way, this is impractical 
since the coefficients quickly become extremely compli- 
cated. Hence we have computed the recursion coefficients 
numerically. When doing so, one encounters another 
problem. The Stieltjes algorithm is extremely-ill con- 
ditioned [IHl , implying that small errors blow up quickly. 
To avoid this, we followed Ref. [isj, and performed the 
Stieltjes algorithm using arbitrary precision arithmetic. 
This can be done using for example the computer pro- 
gram Mathematica. In order to compute all recursion 
coefficients up to = 128 with 22 digits accuracy (so 
that it fits in double precision) we had to use 10,000 dig- 
its precision in the Stieltjes procedure. 

In this way we have computed the weights and nodes 
of the Gauss-Maxwell quadrature with p — 1 up to 
N = 128. For iV = 2, 4, 8 and 16 we could compare with 
tables presented in Ref. [i^. We find excellent agree- 
ment, almost up to machine precision accuracy. 



Appendix D: Computation of H 

Here we will compute H, which is defined in Eq. 
The integration over the parts of Hm{^) that are inde- 
pendent of p is straightforward so that we can write 



hQmSji H —Sji, 

' 2M ^' 



where the matrices A and B read 



(Dl) 



(Ah 



A2 



p2 _|_ jy2A2 



(D2) 
(D3) 



First we will compute the matrix {A)ij. We find 



[Ah 



dxw^''^{x)li{x) 

'— + — - 

dx"^ Ax"^ 

dx w{x)li{x) 



w^'\x)l,{x) 



(D4) 



dx^ 



1 \ d ■ 

- - 2a; 2 

X I dx 



lj{x). 



To obtain the last line we have used that w{x) 
X ex-p{—x^). Eq. (|D4p can be rewritten as 



/•OO 

{A)ij = / dxw{x)li{x) 
Jo 



^ _ 2x— - 2 
da;^ da; 



lj{x) 



+ / dx w{x)li{x) — 



-;,(a;)-Z;(0) 



+ l'j{0) dxw{x)-[k{x) -k{0)] 
Jo ^ 

/•OO 

+ k{0)l'j{0) / da; exp(-x2). 
Jo 



(D5) 



This form has the advantage that the integrands of the 
first three terms are products of the weight function w(x) 
and a polynomial of degree less than 2N. Hence we can 
evaluate these terms exactly using the Gauss-Maxwell 
quadrature. The last term of Eq. (|D5I) can also be com- 
puted analytically and is equal to li{0)lj{0)y/n/2. 

To proceed we will first evaluate the first and second 
order derivatives of Ij (x) at the nodes Xi which become 
(see also Ref. (H) 



1 



64 



1 C,: 



1 



(D6) 



'^1 ^ ' 



■(1-4), 



dx^ 



where 



wl Cj 



2^. 



X'l Xj i^Xi •^j)'^ 



E 



N 



Q = ^/w~ 



(1-4), (D7) 

(D8) 
(D9) 
(DIO) 
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Furthermore li (0) and E (0) are given by 



kiO) 



1 ^ 

^ n 



^ 1 
= ^- 

n—l.n^j 

Putting everything together we find that 
1 
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2a;, j - 2 
2 



■ TV 

E 



1 



/,(0)/,(0). 



(Dll) 
(D12) 
(D13) 



(1 - 5^,) 



(D14) 



times 'w{x). We can then easily integrate these terms us- 
ing the Gauss-Maxweh quadrature. We can rewrite (-B) y 
in the hmit 77 — >■ as 

= / Axw{x) [k{x) - [l,{x) - /,(0)] \ 

JQ ^ 

poo , 

+ k (0) / dx w{x) [I J (x) - (0) - xl'j (0)] — 
Jo ^ 

r°° 1 

+ I, (0) / dx w{x) [k {x) - k (0) -xl[{Q)] — 

JQ 2; 

+ [;.(0);;(0) + ^:(0);,(0)] / da: exp(-a:2) 

Jo 



+ li{0)ljiO) / dx 



(D16) 



The first three terms of the last equation can be com- 
puted using the Gauss-Maxwell quadrature. The next 
term can be evaluated analytically. The last integral can 
be computed analytically for small t]. We obtain 



The matrix A is symmetric as follows from Eq. (|D4I) . 
although this is not directly clear from the last equation. 

Now let us compute the matrix {B)ij. Expressed in 
terms of an integral over the basis functions this matrix 
reads 



.00 ^ 
{B)ij = / dxw{x)k{x)lj{x)^ — 
Jo X + 



(D15) 



We are interested in (-B)y in the limit of small ry. There- 
fore in the rest of the calculations, we will drop terms 
that are of order 77 and higher. Like in the calculation 
for {A)ij we rewrite the integrand in such a way that we 
get terms which are a polynomial of degree less than 2N 



1 ^ 

1 V ^ Wn 

^lE + log(ry) + 2^ 



;,(0);,(0), (D17) 



where 7^ denotes the Euler-Mascheroni constant. 

The last equation shows that has a logarithmic 

singularity for 77 = 0. For that reason we had to regular- 
ize the centrifugal potential in Eq. (|55p . An alternative 
way of treating such singularity in the DVR method is 
discussed in Ref. [57 1. 
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